{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merging radar rain images and wind predictions in a deep learning model applied to rain nowcasting : Training procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import logging\n",
    "import argparse\n",
    "from tqdm import tqdm\n",
    "\n",
    "import torch\n",
    "import numpy as np\n",
    "import torch.nn as nn\n",
    "from torch import optim\n",
    "from torch.utils.data import DataLoader\n",
    "from torch.utils.tensorboard import SummaryWriter\n",
    "\n",
    "import samplers\n",
    "import utilities as utils\n",
    "from unet_model import UNet\n",
    "from dataset import MeteoNetDataset\n",
    "from eval import eval_net_and_persistance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Path to data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directory where rain data is stored\n",
    "rain_dir = r'data/Rain'                      # <- TO COMPLETE\n",
    "# Directory where wind data are stored\n",
    "U_dir = r'data/U'                            # <- TO COMPLETE\n",
    "V_dir = r'data/V'                            # <- TO COMPLETE\n",
    "# Path to rain radar coordonates file\n",
    "coord_file = r'data/radar_coords_NW.npz'          # <- TO COMPLETE\n",
    "# Boolean to save or not the model at the end of each epoch\n",
    "save_cp = False\n",
    "# Directory to save checkpoints if save_cp==True\n",
    "ckp_dir = r''\n",
    "#Architecture constants\n",
    "Matrix_path = os.path.join(\"Matrix\", \"Matrix\")  # Relative path to matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of time steps to use for inputs \n",
    "temporal_length_inputs = 12 #1h input\n",
    "# Number of time steps from the first input to the target\n",
    "temporal_length = 18 # Prediction at 30m. Set 24 for 1h.\n",
    "# Thresholds in mm/h and mm\n",
    "thresholds_in_mmh = [0, 0.1, 1, 2.5]  # CRF over 1h\n",
    "thresholds_in_cent_mm = [100*k/12 for k in thresholds_in_mmh] # CRF over 5 minutes in 1/100 of millimeters like MN data\n",
    "# Input and output channels\n",
    "n_channels = 3*temporal_length_inputs # Factor 3 for rain + U + V\n",
    "n_classes = len(thresholds_in_cent_mm)-1\n",
    "# Proportionnal to the number of weights in the model\n",
    "model_Size = 1 # 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters of the training procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 32\n",
    "epochs = 20\n",
    "percentage_sampling = None # see paper section oversampling\n",
    "# Learning rate\n",
    "lr = 0.0008\n",
    "# Cuda, number of loader worker processes.\n",
    "num_workers = 0\n",
    "# Weight decay for L2 regularization\n",
    "wd = 1e-5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(net,\n",
    "          rain_dir,\n",
    "          U_dir,\n",
    "          V_dir,\n",
    "          ckp_dir,\n",
    "          device,\n",
    "          epochs=20,\n",
    "          batch_size=256,\n",
    "          lr=0.0008,\n",
    "          save_cp=True,\n",
    "          num_workers=0,\n",
    "          n_val_round_by_epoch=3,\n",
    "          wd=1e-5,\n",
    "          percentage_sampling=0.8):\n",
    "    \n",
    "    train = MeteoNetDataset(rain_dir=os.path.join(rain_dir,'train'), \n",
    "                            U_dir=os.path.join(U_dir,'train'), \n",
    "                            V_dir=os.path.join(V_dir,'train'),\n",
    "                            temporal_length_inputs=temporal_length_inputs, \n",
    "                            temporal_length=temporal_length, \n",
    "                            temporal_stride=temporal_length_inputs, \n",
    "                            thresholds=thresholds_in_cent_mm, \n",
    "                            Matrix_path=Matrix_path)\n",
    "    val = MeteoNetDataset(rain_dir=os.path.join(rain_dir,'val'), \n",
    "                          U_dir=os.path.join(U_dir,'val'), \n",
    "                          V_dir=os.path.join(V_dir,'val'),\n",
    "                          temporal_length_inputs=temporal_length_inputs, \n",
    "                          temporal_length=temporal_length, \n",
    "                          temporal_stride=temporal_length_inputs, \n",
    "                          thresholds=thresholds_in_cent_mm, \n",
    "                          Matrix_path=Matrix_path)\n",
    "    thresholds_normalized = np.log(np.array(thresholds_in_cent_mm)+1)/train.norm_factor\n",
    "    \n",
    "    train_sampler, real_percentage = samplers.oversample_xpercent_rainy_tiles(train, p=percentage_sampling, \n",
    "                                                                              above_class=len(thresholds_in_cent_mm)-1)\n",
    "    train_loader = DataLoader(train, batch_size=batch_size, sampler=train_sampler, num_workers=num_workers, pin_memory=True)\n",
    "    # Custom sampler to consider only defined data\n",
    "    val_sampler = samplers.CustomSampler(samplers.indices_except_undefined_sampler(val))\n",
    "    val_loader = DataLoader(val, batch_size=batch_size, sampler=val_sampler, num_workers=num_workers, pin_memory=True, drop_last=True)\n",
    "    \n",
    "    writer = SummaryWriter(comment=f'-LR_{lr}_BS_{batch_size}_E_{epochs}_T{temporal_length}_WD_{wd}_p{percentage_sampling}')\n",
    "    writer.add_scalar('Number_of_parameters', net.get_nb_params())\n",
    "    info = f'''Starting training:\n",
    "        Epochs:                {epochs}\n",
    "        Learning rate:         {lr}\n",
    "        Batch size:            {batch_size}\n",
    "        Weight decay:          {wd}\n",
    "        Number batch train :   {len(train_loader)}\n",
    "        Number batch val :     {len(val_loader)}\n",
    "        Training size:         {len(train)}\n",
    "        Validation size:       {len(val)}\n",
    "        Checkpoints:           {save_cp}\n",
    "        Device:                {device.type}\n",
    "        Number of parameters:  {net.get_nb_params()}\n",
    "        Thresholds(mm/h):      {thresholds_in_mmh}\n",
    "        Temporal length:       {temporal_length}\n",
    "        Normalization factor:  {train.norm_factor}\n",
    "        Percentage sampling:   {percentage_sampling}\n",
    "        Real percantage :      {real_percentage}\n",
    "    '''\n",
    "    logging.info(info)\n",
    "    writer.add_text('Description', info)\n",
    "\n",
    "    global_step = 0\n",
    "    optimizer = optim.Adam(net.parameters(), lr=lr, weight_decay=wd)\n",
    "    criterion = nn.BCEWithLogitsLoss()\n",
    "\n",
    "    for epoch in range(epochs):\n",
    "        net.train()\n",
    "        epoch_loss = 0\n",
    "        n_batch = len(train_loader)\n",
    "        \n",
    "        with tqdm(total=n_batch, desc=f'Epoch {epoch + 1}/{epochs}', unit='batch') as pbar:\n",
    "            for batch in train_loader:\n",
    "                imgs = batch['inputs']  #BC(temp)HW\n",
    "                true_imgs = batch['target']   #BClsHW\n",
    "                assert imgs.shape[1] == net.n_channels, \\\n",
    "                    f'Network has been defined with {net.n_channels} input channels, ' \\\n",
    "                    f'but loaded images have {imgs.shape[1]} channels. Please check that ' \\\n",
    "                    'the images are loaded correctly.'\n",
    "                   \n",
    "                imgs = imgs.to(device=device, dtype=torch.float32)\n",
    "                true_imgs = true_imgs.to(device=device, dtype=torch.float32)\n",
    "                \n",
    "                imgs_pred = net(imgs)   #BClsHW\n",
    "                loss = criterion(imgs_pred, true_imgs)\n",
    "                epoch_loss += loss.item()\n",
    "                pbar.set_postfix(**{'loss (epoch)': epoch_loss})\n",
    "                optimizer.zero_grad()\n",
    "                loss.backward()\n",
    "                nn.utils.clip_grad_value_(net.parameters(), 0.1)\n",
    "                optimizer.step()\n",
    "                pbar.update()\n",
    "                \n",
    "                if ((global_step % (n_batch//n_val_round_by_epoch))==0):\n",
    "                    Scores = eval_net_and_persistance(net=net, thresholds=thresholds_normalized, loader=val_loader, device=device)\n",
    "                    # Log scores\n",
    "                    for key in Scores.keys():\n",
    "                        score = Scores[key]\n",
    "                        if type(score)==type(dict()):\n",
    "                            writer.add_scalars(f'{key}_val', score, global_step)\n",
    "                        else :\n",
    "                            writer.add_scalar(f'{key}_val', score, global_step)\n",
    "                            \n",
    "                    # Save images of batch inputs//target, prediction, persistance, (pred-target) and (pers-target) on the last epoch\n",
    "                    if (epoch==epochs-1):    \n",
    "                        # we select rain, U and V channels (BCHW)\n",
    "                        utils.writer_add_batch_rain(rain_channels=imgs[:,:12,:,:],text=\"Rain_inputs\",writer=writer)\n",
    "                        utils.writer_add_batch_wind(wind_channels=imgs[:,12:24,:,:],text=\"U_inputs\",writer=writer,isU=True)\n",
    "                        utils.writer_add_batch_wind(wind_channels=imgs[:,24:,:,:],text=\"V_inputs\",writer=writer,isU=False)\n",
    "                        utils.writer_add_comparison(imgs,true_imgs,imgs_pred,\n",
    "                                                    text=\"Target_Prediction_Persistance_Prediction-Target_Persistance-Target\",\n",
    "                                                    writer=writer,thresholds=thresholds_normalized, temporal_length_inputs=temporal_length_inputs)\n",
    "        \n",
    "                global_step += 1\n",
    "        losses = {\"lossTot\" : epoch_loss}\n",
    "        writer.add_scalars('Loss_train', losses, epoch)\n",
    "        \n",
    "        if save_cp:\n",
    "            try:\n",
    "                os.mkdir(ckp_dir)\n",
    "                logging.info('Created checkpoint directory')\n",
    "            except OSError:\n",
    "                pass\n",
    "            torch.save(net.state_dict(),\n",
    "                       ckp_dir + f'CP_epoch{epoch + 1}.pth')\n",
    "            logging.info(f'Checkpoint {epoch + 1} saved !')\n",
    "    writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s')\n",
    "    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    logging.info(f'Using device {device}')\n",
    "    \n",
    "    net = UNet(n_channels, n_classes, bilinear=True, n=model_Size)\n",
    "    net.to(device=device)\n",
    "    \n",
    "    try:\n",
    "        train(net=net,\n",
    "              rain_dir=rain_dir,\n",
    "              U_dir=U_dir,\n",
    "              V_dir=V_dir,\n",
    "              ckp_dir=ckp_dir,\n",
    "              epochs=epochs,\n",
    "              batch_size=batch_size,\n",
    "              lr=lr,\n",
    "              device=device,\n",
    "              num_workers=num_workers,\n",
    "              wd = wd,\n",
    "              percentage_sampling=percentage_sampling,\n",
    "              save_cp=save_cp)\n",
    "        \n",
    "    except KeyboardInterrupt:\n",
    "        try:\n",
    "            sys.exit(0)\n",
    "        except SystemExit:\n",
    "            os._exit(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The tensorboard extension is already loaded. To reload it, use:\n",
      "  %reload_ext tensorboard\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "ERROR: Timed out waiting for TensorBoard to start. It may still be running as pid 6188."
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%load_ext tensorboard\n",
    "logs_dir = './runs'\n",
    "%tensorboard --logdir {logs_dir}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An example of prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Using device cuda\n",
      "INFO: Model loaded from D:\\Repositories\\rain-nowcasting-with-fusion-of-rainfall-and-wind-data-article\\weigths\\model30m.pth\n",
      "INFO: Creating dataset\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAIlCAYAAABRrIC2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdf7xddX3n+9eHUEEpWFJIJkIQKtHeUCvXk4u1gFfrVWLHGOpUwZYpbblNx0JL6S/Dve3IOLbj7bRmmA5a08okHaoEZ+RBGCmUZqqAd/oAjmMVf6TkAkKaSKCxgo6Cwc/946ydrOzsH+fss/dea6/9ej4e57H3+u691/ks4LxZn/1dPyIzkSRJkiRV66iqC5AkSZIk2ZxJkiRJUi3YnEmSJElSDdicSZIkSVIN2JxJkiRJUg3YnEmSJElSDRxddQEazEknnZSnn346s7Ozi17XzMxM19eGsf4RejIzT666CEmHRIT3Z5ljPkk1MspsmpmZqfv+UpnZVHPhfc4m05o1a3KYQdD+30FEDG3dIzSbmWuqLkLSITZnB5lPUo2YTQeZTTXnzNmEGlVjNiFNmSRJktQ4nnM25WzMJEmSpHpw5mzKzachazVwNm+SJEnS6Dhzpnnz/ERJkiRpdJw507w4ayZJkiSNljNn6svGTJIkSRo9mzNJkiRJqgEPa5ximUlEdDyXzNkySXXhRYkkSdPC5mzKdWvM3BmSVLVO+dT6UkmSquKX2holD2scgYi4PiL2RcQDpbGlEXFnRDxYPJ5Yeu3qiNgVETsj4oIx1nnwp8zGTGqmScmm0u8/LIe8YqzUTJOaTd32o6TFsDkbjS3A2raxjcCOzFwF7CiWiYjVwMXAWcVnPhARS8ZX6pxyuBg0UmNtYcKyqajFnSCp2bYwgdkkjYLN2Qhk5l3A/rbh9cDW4vlW4MLS+I2Z+UxmPgzsAs4ZS6EduOMjNdckZ5Ok5jKbpENszsZneWbuBSgelxXjpwCPld63uxg7QkRsiIj7I+L+kVYqaZqYTZLqyGzSVPKCINXrNFXV8cSKzNwMbAaICE++kDRKZpOkOjKb1GjOnI3P4xGxAqB43FeM7wZWlt53KrBnzLVJml5mk6Q6Mps0lWzOxmc7cGnx/FLgltL4xRFxTEScAawC7q2gPknTyWySVEdmk6aSzVkXEfH/zGesy2c/Cvx34GURsTsiLgPeB7whIh4E3lAsk5lfAG4CvgjcDlyemc8NZysk6RCzSVIdmU3SIeF9YzqLiM9k5ivbxj6XmT9cVU1lEZHejJXZzFxTdRGSDvG8joPMJ6lGzKaDzKaac+asTUS8MyI+z9y3N58r/TwMfK7q+lpmZmYAb8oqqV5mZmbMJUmSBuTVGo/0EeAvgH9DccPDwtOZ2X4PDklSyezsbNUlSJI0sZw5a5OZX8/MRzLzHcxdDejHMvMrwFHFiae1MDs7S0RM+2GNkmrIXJIkaTA2Z11ExLuBdwFXF0PPA26oriJJkiRJTWZz1t1PAG8BvgmQmXuA4yutSJIkSVJj2Zx192zOndWeABFxXMX1SJIkSWowm7PuboqIDwHfFxG/APwV8CcV1yRJkiSpobxaYxeZ+QcR8QbgKeBlwL/MzDsrLkuSJElSQ9mc9VA0YzZkkiRJkkbO5qyLiHia4nyzkq8D9wO/npkPjb8qSZIkSU1lc9bd+4E9zN2UOoCLgX8C7ASuB15bWWWSJEmSGscLgnS3NjM/lJlPZ+ZTmbkZ+PHM3AacWHVxkiRJkprF5qy770bE2yPiqOLn7aXX2g93lCRJkqRFsTnr7qeBfw7sAx4vnl8SEc8HrqiyMEmSJEnN4zlnHUTEEuCdmbmuy1vuGWc9kiRJkprPmbMOMvM5YKbqOiRJkiRND2fOuvsfEbEd+BjwzdZgZn58MSuNiEeAp4HngAOZuSYilgLbgNOBR4C3Z+bXFvN7JGkhzCZJdWU+aZo4c9bdUuAfgB8D1hU/bx7Sul+XmWdn5ppieSOwIzNXATuKZUkaN7NJUl2ZT5oKzpx1kZk/N8Zft55D903bCnwSeNcYf78kdWI2Saor80mNZHPWRUQcC1wGnAUc2xrPzJ9f5KoT+MuISOBDxf3Tlmfm3mL9eyNi2SJ/x8KLykN3B4iIcf96SdWrZTZJEuaTpojNWXf/CfgycAHwHuYurf+lIaz33MzcU4TInRHx5fl+MCI2ABuGUEOndY9itSPRaiQnqWZpAtQymySJAfPJbNIk8pyzNhHRaljPzMzfAb6ZmVuBfwq8fLHrz8w9xeM+4GbgHODxiFhR/P4VzN1brdNnN2fmmtLx1lMjMw/+SBo+s2lw5Xwyp6ThGzSfpj2bWsymyWJzdqR7i8fvFI//GBE/BLyQuSsCDSwijouI41vPgTcCDwDbgUuLt10K3LKY3zPpyuFhmEijZzYNznySRst8GowN2eTysMbuNkfEicBvMxcA3wv8ziLXuRy4uTgc72jgI5l5e0TcB9wUEZcBjwJvW+TvmSidgsMwkcbKbOrAHJJqwXzqIjOJiHlnlaeDTIbwfz6Hi4jdwPvbh4vHzMz21ypRnBRbW53+u2qFwrD+m4uI2Wk/VEGqm7pnU0trp6bba4MqrdN8kmpkUrJpIfplVZeMM5tqzpmzIy1hbpas03/RjfvDHoVuYTGMpsxvfSQtViuLzCRJddL+pdGgGWUuTTabsyPtzcz3VF3EpBjXJfgNGknDMMyjRcwlScO0kEMU1Vw2Z0fy/7YlCwmJYQeKOz6SFqrXF0bOlEmqsyEdUq0JZ3N2pNdXXUDVqv7WxoCR1E0VXxiZSZLmq1vujOC8+6GsR/Vjc9YmM/dXXcMoVd14lRkskrqpKqvMJUkLMd8rJs734h3l885anzGXpovNWYPVoREzUCR1YwMmaVJ1aqIWuz44PJ/MqulkcyZJqkSvHY9RNm69LqMvSf0s9oqvnWbazCS12Jw1QK/jm4f9zc58GTKSyubTEI3igh3uAEmqWqfc6XQYowRwVNUFaDAzMzNkZs+dmfJx0Ivd6YmIrj+d3iNJZcPIhV451C17+r0uSfPVaX9qkFxq/7xUZnM2oWZnZ/u+Z5h/8L2aO3d6JI1KK1vMGElV6/SFtDRsNmcTamZm5uDz9nDoNKO1WAaQpHGzMZNUJ60jkswkjZLnnDXAYo5X7vY5L98qaVw6zcybPZLqxlzSONicTbh+Vzsrf/Ncbrj6NXQGkKRxMW8kSZpjc9ZgnQ537PaaJEmSpGp5ztmEmp2dtcGSJEmSGsTmTJIkSZJqwOZMkiRJkmrA5qwmImJtROyMiF0RsbHqeiQJzCZJ9WU+qYlszmogIpYA1wFvAlYD74iI1dVWNV6ZOe8fSeNhNkmqK/NJTWVzVg/nALsy86HMfBa4EVhfcU1jsZCGyxs/SmM3tdkkqfbMJzWSl9Kvh1OAx0rLu4FXtb8pIjYAG4rFZ4AHRl/aaBXN1knAkwN8/MXDrUZSm6nNpsKg2QTmkzRqffOpwdkE7js1ls1ZPXSaDjpiOikzNwObASLi/sxcM+rCxqFJ2yI1jNnUkG2RGqhvPjU1m6B526NDPKyxHnYDK0vLpwJ7KqpFklrMJkl1ZT6pkWzO6uE+YFVEnBERzwMuBrZXXJMkmU2S6sp8UiN5WGMNZOaBiLgCuANYAlyfmV/o87HNo69sbJq0LVJjmE2N2hapUQbIp6b9PTdte1QIL00uSZIkSdXzsEZJkiRJqgGbM0mSJEmqAZuzCRMRayNiZ0TsioiNVdfTT0RcHxH7IuKB0tjSiLgzIh4sHk8svXZ1sW07I+KCaqqWNAjzSVIdmU2aJDZnEyQilgDXAW8CVgPviIjV1VbV1xZgbdvYRmBHZq4CdhTLFNtyMXBW8ZkPFNssqebMJ0l1ZDZp0ticTZZzgF2Z+VBmPgvcCKyvuKaeMvMuYH/b8Hpga/F8K3BhafzGzHwmMx8GdjG3zZLqz3ySVEdmkyaKzdlkOQV4rLS8uxibNMszcy9A8bisGG/K9knTqCl/v+aT1CxN+ds1m6aEzdlkiQ5jTboXQtO3T2qypv/9Nn37pKZq+t9u07dv6ticTZbdwMrS8qnAnopqWYzHI2IFQPG4rxhvyvZJ06gpf7/mk9QsTfnbNZumhDehniARcTTwd8DrgYdG8TtmZmYOPp+dnR3FrximJzPz5KqLkHRYPp1RdS01YT5JNWA2HcFsqrmjqy5A85eZByLiCuCOEa2/43hEpxnzWvhK1QVImlPKp09UXUtNmE9SDZhNRzCbas7DGidMZt6WmS8d8+8c56+TNKEy87aqa5CkdmaTJonNmQ6q8QyZJEmS1Hg2Z5oXZ88kSZKk0bI5kyRJkqQasDmTJEmSpBqwOdNBvQ5d9Hw0SZIkabRsztSXjZkkSZI0ejZn6snGTJIkSRoPm7MRiIjrI2JfRDxQGlsaEXdGxIPF44ml166OiF0RsTMiLqim6iMbMRszqVkmNZskNZvZJB1ic9ZHRJxbhMLfRcRDEfFwRDzU52NbgLVtYxuBHZm5CthRLBMRq4GLgbOKz3wgIpYMdSMWICIO/khqnC1MaDZJarQtmE0SYHM2Hx8G3g+cB/xvwJrisavMvAvY3za8HthaPN8KXFgavzEzn8nMh4FdwDnDKV2SDjGbJNWR2SQdcnTVBUyAr2fmXwxhPcszcy9AZu6NiGXF+CnA35Tet7sYk6RxMJsk1ZHZpKlkc9bfX0fEvwU+DjzTGszMzwxp/Z2OH+x4TfuI2ABsGNLvlaRezCZJdWQ2qdFszvp7VfG4pjSWwI8tcD2PR8SK4tufFcC+Ynw3sLL0vlOBPZ1WkJmbgc0AEdH9pmSSNH9mk6Q6Mps0lTznrI/MfF2Hn4U2ZgDbgUuL55cCt5TGL46IYyLiDGAVcO/iK5ekeTGbJNWR2aSp5MxZHxHxQuDdwGuKoU8B78nMr/f4zEeB1wInRcTu4vPvA26KiMuAR4G3AWTmFyLiJuCLwAHg8sx8bkSbI2mKmU2S6shskg6JTGd5e4mI/wI8wKErBv1z4BWZ+dbqqnJ6vjCbmWv6v03SuJhNB5lPUo2YTQeZTTXnzFl/L8nMf1Za/lcR8dnKqpEkSZLUSJ5z1t+3IuK81kJEnAt8q8J6JEmSJDWQM2f9vRPYWpx7FszdJPFnK61IkiRJUuPYnPWRmZ8FXhERJxTLT1VckiRJkqQGsjnrIiIuycwbIuLX2sYByMz3V1KYJEmSpEayOevuuOLx+A6vecUfSZIkSUNlc9ZFZn6oePpXmfnp8mvFRUEkSZIkaWi8WmN/fzTPMUmSJEkamDNnXUTEq4EfBU5uO+/sBGBJNVVJkiRJaiqbs+6eB3wvc/+MyuedPQX8ZCUVSZIkSWosm7MuMvNTwKciYktmfqXqeiRJkiQ1m+ec9fenEfF9rYWIODEi7qiyIEmSJEnNY3PW30mZ+Y+thcz8GrCswnokSZIkNZDNWX/fjYjTWgsR8WK8z5kkSZKkIfOcs/7+b+CeiPhUsfwaYEOF9UiSJElqIJuzPjLz9oh4JfAjQABXZeaTFZclSZIkqWE8rLGLiPjB4vGVwGnAHuDvgdOKMUmSJEkaGmfOuvt14BeAP+zwWgI/NshKI+IR4GngOeBAZq6JiKXANuB04BHg7cWFRyRpLMwmSXVlPmmaRKbXthinImDWlA+NjIjfB/Zn5vsiYiNwYma+q896/BcHs5m5puoipCYwm4bOfJKGZBj5ZDYdZDbVnDNnXUTEW3u9npkfH+KvWw+8tni+Ffgk0HMHSJLGwGySVFfmkxrJ5qy7dcXjMuBHgf9WLL+OuQAYtDlL4C+Lb3A+lJmbgeWZuRcgM/dGhPdR66I80xsRFVYiNY7ZJKmuzCdNDZuzLjLz5wAi4r8Cq1sBEBErgOsWsepzM3NPESJ3RsSX5/vBiNjAlF7G38NvpZEzmyTV1UD5ZDZpEnm1xv5ObzVmhceBlw66sszcUzzuA24GzgEeL5q+VvO3r8tnN2fmmmk5VjgzD/5IGi2zaTDlnDKzpNEYNJ+mOZs0uWzO+vtkRNwRET8bEZcCnwD+epAVRcRxEXF86znwRuABYDtwafG2S4FbFl/2ZHLnRho/s2lhzClpfMynwfil0eTysMY+MvOKiPgJ4DXF0ObMvHnA1S0Hbi7OlToa+Ehxk+v7gJsi4jLgUeBti617UhgWUi2YTX2YVVJlzKd5MqeawUvpz0NEvBhYlZl/FREvAJZk5tMV11Tbf3Hj+m8qIrwcrFQzdc6m+RhWfplPUr1MejZ1stC8Khpcs6nmnDnrIyJ+gbmTSZcCLwFOAf4YeH2VddWNTb6kOjOjJNWdOSXwnLP5uBw4F3gKIDMfZO7y+lPN45gl1V1VGeVtPiTNx7j2pSLCXJogzpz190xmPtv6jzoijmbufhtTw+ZL0qSoKq/c8ZE0H+PIKPNostmc9fepiPi/gOdHxBuAXwJurbimkal7I2bgSKpDTplFkuZjjOfhj+X3aPQ8rLG/dwFPAJ8HfhG4DfjtSisasjofntiaindKXppudcgps0hSL1Vcvt5Mah5nznqIiKOAz2XmDwF/UnU9w1THRqzFkJHUUnVWmUeS+hlnTplJzWdz1kNmfjci/jYiTsvMR6uuZzGq3sHpxaCRVFZFXplDkqpmDglszuZjBfCFiLgX+GZrMDPfUl1J81OnhszAkdTLqPPKDJI0KuaLhsnmrL9/VXUBk8JwkjQo80OSJJuzriLiWOBfAGcydzGQD2fmgWqrWpiIWNC30e4cSZIkSdWxOetuK/Ad4G7gTcBq4MpKKxpAueEqN2o2YpIkSVK92Jx1tzozXw4QER8G7q24nkWzIZMkSZLqy/ucdfed1pNJO5xRkiRJ0uRx5qy7V0TEU8XzAJ5fLAeQmXlCdaVJkiRJahqbsy4yc0nVNUiSJEmaHh7WKEmSJEk1YHMmSZIkSTVgc1YTEbE2InZGxK6I2Fh1PZIEZpOk+jKf1EQ2ZzUQEUuA6zh0P7V3RMTqaquSNO3MJkl1ZT6pqWzO6uEcYFdmPpSZzwI3AusrrkmSzCZJdWU+qZG8WmM9nAI8VlreDbyq/U0RsQHYUCw+Azww+tLG4iTgyQE+9+JhFyLpMGbTYNkE5pM0an3zqcHZBO47NZbNWT1Eh7E8YiBzM7AZICLuz8w1oy5sHJq0LVLDmE0N2RapgfrmU1OzCZq3PTrEwxrrYTewsrR8KrCnolokqcVsklRX5pMayeasHu4DVkXEGRHxPOBiYHvFNUmS2SSprswnNZKHNdZAZh6IiCuAO4AlwPWZ+YU+H9s8+srGpknbIjWG2dSobZEaZYB8atrfc9O2R4XIPOL0AUmSJEnSmHlYoyRJkiTVgM2ZJEmSJNWAzdmEiYi1EbEzInZFxMaq6+knIq6PiH0R8UBpbGlE3BkRDxaPJ5Zeu7rYtp0RcUE1VUsahPkkqY7MJk0Sm7MJEhFLgOuANwGrgXdExOpqq+prC7C2bWwjsCMzVwE7imWKbbkYOKv4zAeKbZZUc+aTpDoymzRpbM4myznArsx8KDOfBW4E1ldcU0+ZeRewv214PbC1eL4VuLA0fmNmPpOZDwO7mNtmSfVnPkmqI7NJE8XmbLKcAjxWWt5djE2a5Zm5F6B4XFaMN2X7pGnUlL9f80lqlqb87ZpNU8LmbLJEh7Em3Quh6dsnNVnT/36bvn1SUzX9b7fp2zd1vM/ZBImIVwPXZOYFEeG/OHgyM0+uughJh/IJeGPFpdSF+STVgNl0BLOp5pw5myz3Aasi4oyZmRkykylvrr9SdQGSDroPWAWYTXPMJ6keDmaTALOp9mzOJkhmHgCuAO6YnZ0lIogId4IkVa6UTwezSZKqVs4maRLYnE2YzLwtM1/aPnNmgyapapl5WzmbnEGTVAedskmqK5uzBjFsJNWNM2iSJM3f0VUXoOFo7QBlpjtDkipnDkmqK/NJdWZzNuHaA8bAkSRJOlzrXH2p7jyssY+IeEFE/E5E/EmxvCoi3lx1XZJUV+4ASZI0GJuz/v4j8Azw6mJ5N/De6so5xB0gSXUzOztbdQmSJE0sm7P+XpKZvw98ByAzv0Xnu7GPVacdIK9CJKmuzCZJkvrznLP+no2I5wMJEBEvYW4mrTbc4ZFUZ+ULFXnRIkmSunPmrL93A7cDKyPiz4EdwG/1+kBEXB8R+yLigdLY0oi4MyIeLB5PLL12dUTsioidEXHBQgts3fC1/CNJ7cadTaX1dHwuSVBdNkl1ZHPWR2beCbwV+Fngo8CazPxkn49tAda2jW0EdmTmKuYavI0AEbEauBg4q/jMByJiyZDKl6SyLZhNkupnC2aTBNic9RURPwEcyMxPZOZ/BQ5ExIW9PpOZdwH724bXA1uL51uBC0vjN2bmM5n5MLALOGdoGyBJBbNJUh2ZTdIhNmf9vTszv95ayMx/ZO5Qx4Vanpl7i3XsBZYV46cAj5Xet7sY62lmZuawC4B43pmkAZlNkupoqNkkTQqbs/46/TMa5oVUOp2A0XFvJiI2RMT9EXH/E088ceSH3BGSNDxDyyZJGqKBsmnENUlD49Ua+7s/It4PXMfcH/8vA4PcyOfxiFiRmXsjYgWwrxjfDawsve9UYE+nFWTmZmAzQETYhUkahqFnkxf9kDQE7jdpKjlz1t8vA88C24CPAd8GLh9gPduBS4vnlwK3lMYvjohjIuIMYBVw73xW2H4FNK/UKGkAQ88mSRoCs0lTyZmzPjLzmxRXCJqviPgo8FrgpIjYzdw5au8DboqIy4BHgbcV6/9CRNwEfBE4AFyemc8t4HctpDRJU2yc2SRJ82U2SYeE5yj1FhEvBX4DOJ1SM5uZP1ZVTeD0fGE2M9dUXYSkQ8ymg8wnqUbMpoPMpppz5qy/jwF/DPwp4DczkiRJkkbC5qy/A5n5waqLkCRJktRsXhCkv1sj4pciYkVELG39VF3UtPNwXEmSJDWNM2f9ta4U9JulsQR+oIJaxFxj5oVQJEmS1DQ2Z31k5hlV16DD2ZhJkiSpiWzO5iEifghYDRzbGsvMP6uuoungDJkkSZKmic1ZHxHxbubuvbEauA14E3APYHM2BjZokiRJmhZeEKS/nwReD3w1M38OeAVwTLUlNVdmHvyRJEmSponNWX/fyszvAgci4gRgH14MRJIkSdKQeVhjf/dHxPcBfwLMAt8A7q22pOZpzZR5CKMkSZKmlc1ZH5n5S8XTP46I24ETMvNzVdbURDZlkiRJmnYe1thHROxoPc/MRzLzc+UxSZIkSRoGZ866iIhjgRcAJ0XEiUBraucE4EWVFSZJkiSpkWzOuvtF4FeZa8RmOdScPQVcV1VRkiRJkprJ5qyLzLwWuDYifjkz/6jqeiRJkiQ1m+ec9ffViDgeICJ+OyI+HhGvrLooSZIkSc1ic9bf72Tm0xFxHnABsBX4YMU1TSxvMi1JkiR1ZnPW33PF4z8FPpiZtwDPG3RlEfFIRHw+Ij4bEfcXY0sj4s6IeLB4PHEIdddO+V5mXjpfqpdpziZJ9WY+aZrYnPX39xHxIeDtwG0RcQyL/+f2usw8OzPXFMsbgR2ZuQrYUSw3Uqspc/ZMqqWpzSZJtWc+aSrYnPX3duAOYG1m/iOwFPjNIf+O9cwdLknxeOGQ11+Z9sMYbcqkidLYbJI08cwnNZLNWRcRcULx9Fjgk8A/RMRS4Bng/kWsOoG/jIjZiNhQjC3PzL0AxeOyLjVtiIj7W1P6k85DG6VaMZsk1dVA+WQ2aRJ5Kf3uPgK8mbl7nCWH7nNGsfwDA6733MzcExHLgDsj4svz/WBmbgY2A0RE7aefOs2Q2ZBJtTU12SRp4gyUT2aTJpHNWReZ+ebi8Ywhr3dP8bgvIm4GzgEej4gVmbk3IlYA+4b5O+siIroe0mjTJlVrmrNJUr2ZT5omHtbYR0ScGxHHFc8viYj3R8RpA67ruNI9044D3gg8AGwHLi3edilwy+Irr1anJsxzzaR6mqZs6sfbfUj1Yj5p2jhz1t8HgVdExCuA3wI+DPwn4H8fYF3LgZuLWaKjgY9k5u0RcR9wU0RcBjwKvG0oldecs2VSbUx9NpWbMbNJqhXzKdNcmiLht4O9RcRnMvOVEfEvgb/PzA+3xiquq7b/4ubz39SQQma2dEldSTVQ52xq155VQ975MZ+kGpnkbIKh5pPZVHPOnPX3dERcDVwCvCYilgDfU3FNtTHGRkySFsWLFEmqs1ZGmUvTzXPO+ruIucvnX5aZXwVOAf5ttSXVg42ZpEnQ6RyyiDCfJFWqfI6rjZlanDnroZgluyEz/4/WWGY+CvxZdVXVg42ZpEngbJmkumtdzdpsEtic9ZSZz0XE/4yIF2bm16uup0qemyipzsrfOtuQSaqzTrNkZpRabM76+zbw+Yi4E/hmazAzf6W6kgbTa8q81zc2gzRmhoykKoz4Ih+StGDmkhbC5qy/TxQ/E689DNrDYiFNmMEiadT6zYZ1YjZJqkJ7XrXnltmk+bI56yMzt0bE84HTMnNn1fXMR/ss2LAPSTRgJI1Ct6yyMZNUtfK+Va9Mar3mBT40KK/W2EdErAM+C9xeLJ8dEdurrepI7Vf7aV8eBq9uJmkQnfKp22uDMJskjUqnfat+WplkLmkQzpz1dw1wDvBJgMz8bEScUWVBLaO8SIeBImmxOmVUpwZtvswlSYu12C+CPH9Mo+bMWX8HOlypsfJLF87MzAxtXeVvePymR9JizMzMDGXW3lySVEfmkkbN5qy/ByLip4AlEbEqIv4I+H+rLmoYDBZJdWQuSRq21hdH/XTbNzKXNC42Z/39MnAW8AzwEeDrwK9WWtGQeO8ySePW6xtnvzCSNGrt9xZrLbfnjzNkqornnHUREccC/wI4E/g88OrMPFBtVd0ZHJLqYr55ZG5JGpfZ2dnDGrEys0h14sxZd1uBNcw1Zm8C/qDaco40Dd/odLqqm6T6mp2drboESZImljNn3a3OzJcDRMSHgXsrrucITb2HRr+bYzd1uyVNFrNIkjRsNmfdfaf1JDMP1PV/vnWtaxDOikmqs/ncJLtJmSxJGj8Pa+zuFRHxVPHzNPDDrecR8dSwf1lErI2InRGxKyI29nt/+djpSTaKm2VLGp6FZlOTdbq8v0DlQRUAACAASURBVBcNkKpjPqmJnDnrIjOXjOt3RcQS4DrgDcBu4L6I2J6ZXxxXDVXo1JCVb/DYvqNjAyeN17Rmk6T6M5/UVM6c1cM5wK7MfCgznwVuBNZXXNPItV++ttvlbCVVZiqzSdJEMJ/USDZn9XAK8FhpeXcxNtHmc7jiQpswDx+SxqqR2SSpEcwnNZKHNdZDp07jiK4mIjYAG4rFZ4AHRlnUYi2ggToJeHKA9b14gSVJWphGZtMCzCubujCfpNHqm08NziYYPJ/MppqzOauH3cDK0vKpwJ72N2XmZmAzQETcn5lrxlPeaDVpW6SGMZsasi1SA/XNp6ZmEzRve3SIhzXWw33Aqog4IyKeB1wMbK+4JkkymyTVlfmkRnLmrAaK+6hdAdwBLAGuz8wvVFyWpClnNkmqK/NJTWVzVhOZeRtw2wI+snlUtVSgSdsiNYrZJKmuFphPTft7btr2qBDeO0qSJEmSquc5Z5IkSZJUAzZnEyYi1kbEzojYFREbq66nn4i4PiL2RcQDpbGlEXFnRDxYPJ5Yeu3qYtt2RsQF1VQtaRDmk6Q6Mps0SWzOJkhELAGuA94ErAbeERGrq62qry3A2raxjcCOzFwF7CiWKbblYuCs4jMfKLZZUs2ZT5LqyGzSpLE5myznALsy86HMfBa4EVhfcU09ZeZdwP624fXA1uL5VuDC0viNmflMZj4M7GJumyXVn/kkqY7MJk0Um7PJcgrwWGl5dzE2aZZn5l6A4nFZMd6U7ZOmUVP+fs0nqVma8rdrNk0Jm7PJEh3GmnS5zaZvn9RkTf/7bfr2SU3V9L/dpm/f1LE5myy7gZWl5VOBPRXVshiPR8QKgOJxXzHelO2TplFT/n7NJ6lZmvK3azZNCe9zNkEi4mjg74DXv/CYYx560XHH8djSxd9HfOX+Aweff+mYQ4c4/y/PLB3K+kflG7v2PZmZJ1ddh6RD+XT0CceeceyyE6oup3Lmk1QP48imlfsP1Hp/qcxsqr/J+C9JAGTmgYi4ArjjRccdxw0//uNcedFw/r6u3fYEM2fecNjYF6/Zz/lbrxrK+kfhnnWbvlJ1DZLmtPJp1dEv+MT3b/rpqsupnPkk1UMrm45ddsInzh5hNn3/yNY8XGZT/XlY44TJzNsy86WPLT16aI3ZYeu/5tAPwN3rNg39d0hqpsy8bVK+PZY0PTLztqprkObL5kxcu+2Jnq/boEmSJEmjZ3Omg86b6X4Iow2aJEmSNFoef6K+h0fGNXOPs9ueGMmhlJIkSZJszjRPeQ2cf6uNmSRJkjQqHtaog7odujjXmNX3qo2SJElSE9icSZIkSVIN2Jz1EBFvm8/YpGpdpfHabU9w7bYniGs4+FPmrJmkKrUySpKkpvOcs96uBj42j7GJ1WmHJ6851KD1uoKjJI1Sp3y61gsTSapYp2wylzQsNmcdRMSbgB8HTomIf1966QTgwDw+fz3wZmBfZv5QMbYU2AacDjwCvD0zv1a8djVwGfAc8CuZecfQNqaHcpBcu+0JZnddAsBnLoHZXQaN1DSTkk0trQwqz/JLap5JzSZpFDyssbM9wP3At4HZ0s924IJ5fH4LsLZtbCOwIzNXATuKZSJiNXAxcFbxmQ9ExJLFb8LCXHnRyQfDpvxcUqNsYcKyCQ5lktkkNdYWJjCbpFGwOesgM/82M7cCZ2bm1tLPx1vf2vT5/F3A/rbh9cDW4vlW4MLS+I2Z+UxmPgzsAs4ZzpYsnDs+UnNNcjZJai6zSTrEwxp7+0xEZNvY15mbVXtvZv7DAta1PDP3AmTm3ohYVoyfAvxN6X27i7EjRMQGYAPAMScfv4BfLUldmU2S6shs0lSyOevtL5g7nvkjxfLFQDDXoG0B1g3hd0SHsfaGcG4wczOwGeD4Vcs7vkeShsRsklRHZpMazeast3Mz89zS8ucj4tOZeW5EXLLAdT0eESuKb39WAPuK8d3AytL7TmXunDdJGgezSVIdmU2aSp5z1tv3RsSrWgsRcQ7wvcVi36s2ttkOXFo8vxS4pTR+cUQcExFnAKuAewcvWZIWxGySVEdmk6aSM2e9/Z/A9RHRasieBi6LiOOAf9PtQxHxUeC1wEkRsRt4N/A+4KaIuAx4FHgbQGZ+ISJuAr7IXMN3eWY+N6LtkTTFzCZJdWQ2SYfYnHUQEW8trsx4H/DyiHghEJn5j6W33dTt85n5ji4vvb7L+38X+N2BC5akeTCbJNWR2SQd4mGNnf12eSEzv97WmFVu5f4DXLvtCW/KKqlWWtkkSZIWzuasAdwRklQ35pIkSQvnYY2d/WBEfK7DeACZmT887oK68abRkurksaVHm0uSJA3I5qyzhxnOPcxGxh0gSZIkqVlszjp7NjO/UnURkiRJkqaH55x19un2gYjYXEUhkiRJkqaDzVkHmXlFh+E1Yy9EkiRJ0tSwOZu/fVUXIEmSJKm5bM7mKTPXVl2DJEmSpObygiAdRMStQHZ7PTPfMsZyJEmSJE0Bm7PO/qB4fCvwT4AbiuV3AI9UUZAkSZKkZrM56yAzPwUQEf86M19TeunWiLirorIkSZIkNZjnnPV2ckT8QGshIs4AvPOzJEmSpKFz5qy3q4BPRsRDxfLpwC9WV44kSZKkprI56yEzb4+IVcAPFkNfzsxnqqxJkiRJUjN5WGMPEfEC4DeBKzLzb4HTIuLNFZclSZIkqYFsznr7j8CzwKuL5d3Aexezwoh4JCI+HxGfjYj7i7GlEXFnRDxYPJ64uLIlaWHMJkl1ZT5pmtic9faSzPx94DsAmfktIIaw3tdl5tmZuaZY3gjsyMxVwI5iWZLGzWySVFfmk6aCzVlvz0bE8yluSB0RLwFGcc7ZemBr8XwrcOEIfsfQXbvtiapLkDRaE5lNkqaC+aRGsjnr7d3A7cDKiPhz5r6Z+a1FrjOBv4yI2YjYUIwtz8y9AMXjskX+jrG48iLvKiA1SGOySVLjmE+aGl6tsYuIOAo4EXgr8CPMHc54ZWY+uchVn5uZeyJiGXBnRHx5ATVtADYAHHPy8Yss43CtWbA6N1x3r9t02PIwji+VdFAts0mSGDCfzCZNIpuzLjLzuxFxRWbeBHxiiOvdUzzui4ibgXOAxyNiRWbujYgVwL4un90MbAY4ftXyHFZNMDlNWVxTeuGa9ndKGlRds2lStHIqroHZXZcwU205UqMMmk9mkyaRzVlvd0bEbwDbgG+2BjNz/yAri4jjgKMy8+ni+RuB9wDbgUuB9xWPtyy28EnXtSGTNHRm02DuXrfpYD6Vc2rmzBuqKEdqJPNpca7d9sTBTPKLo8lgc9bbzxePl5fGEviBAde3HLg5ImDun/1Hihtd3wfcFBGXAY8Cbxtw/RPrntm2QxavqaYOaUqZTV1cu+2Jg0cW3DO7ibxmbjyuMaekMTGfuijnU0v5SyOAmTPHW5MWz+ash8w8Y8jrewh4RYfxfwBeP8zfVTft54u1uHMjVW+as6mXe2Y3ze3YzB4a65dZB5u3EdUkTRvzqbfyzBh0z6i8Bs6/9WRwYr/2bM46iIhXMXeM8kuAzwM/n5lfqraq+iufc9HOJkxSHbR/q9yS1ywup1pNmSSNQnsTBt1nxcyjyWZz1tl1wG8AdwFvAf4dcEGlFdVU++GIo9IpaPxmWlI/R1zl9ZrO7xukMeu2A3T+rVdBl6MFJKndfM6zX+zhieffetXiVqCxsTnr7KjMvLN4/rGIuLrSairUfqPpXie6D3t2zG9+JC1U+dvl9gwZRkb1yyV3gCTNR6fTPYY5e28WTS6bs86+LyLe2m05Mz9eQU1jU54Nq+JE0vk0ZX4zLU2n9i+M4PAvjcqZNawvjHplkjtAkso6ZRTM3bKo2xVeB9Epl8yjZrA56+xTwLouywlMbHPWLTSquPTzfGfGDBtp+nRrwkb9hZEzY5IWonXFxHJmddynmu3dkJWvBFt+Xn6tG3OpWWzOOsjMn6u6hmHo1oiVjbopG/TQRINGar75ZBQMP6cGySUzSVJZ+TL2nS7WsVCdGjFn7aeTzVkXEfHBzHxnRFyXmZf3/8RkqPPNUQ0aabq035+nbBg7O92Uv5nuxjyS1E3ri6XW40KzqtPVYXteYEhTxeasg4g4DbgnIrYD2yLitMx8tOq6uul1fHP7TVPHZd7njUmaCp1ultpuGOe79rsIyNy9fq7i/FsHW78ktZvddQnQvUnrtE/UyqJWNppJarE56+x1wGnAy4F7gSXAn1VaUZuV+w/0PSSo1ZiN+gplnabibbwklfVrzGBuB2fmzBsO7uh0Wke3G9p3cv6tV3HeEWPz/rgkHabTflc5286j875Pr9yZTzZquticdZCZWyPiT4FXAe/NzPdWXVO7x5b2/1c3zEvR9zoM6LyZuTByp0fSYrU3Zq0T7Vs7MH7xI6kq5Qt/2FRpVGzOuvswczNn7ykPRsTazLy9mpIOWbn/wMHn7VcJOhgYF83txGTp0q3tun1D3c4pd0nj1soyd4Ik1cF8Ds+WFsvmrIOI+BXgcuBLwNkRcWVm3lK8/HtA5c1ZWb+wOP/Wq5jtcV5at3X2el2ShqXfoUKSVAfmksbB5qyzXwBmMvMbEXE68J8j4vTMvBaISitr0+9qZ+VvnssNV7+GzgCSNC7mjSRJc2zOOluSmd8AyMxHIuK1zDVoL6ZmzVkv7Ts85WV3hiRJkqR6OarqAmrqqxFxdmuhaNTeDJzE3HlotWCDJUmSJDWHzVlnPwN8tTyQmQcy82eA11RT0uHmc7VGSZIkSZPDPfwOMnN3j9c+Pc5aJEmSJE0HZ85qIiLWRsTOiNgVERurrkeSwGySVF/mk5rImbMaiIglwHXAG4DdwH0RsT0zv1htZeNXvopkp8trg+faSeNiNkmqK/NJTWVzVg/nALsy8yGAiLgRWA9MRcC0N2E2ZVJtTHU2Sao180mNZHNWD6cAj5WWdwOvqqiWsfLms1KtTW02Sao980mNZHNWD53unZZHvCliA7ChWHzmnnWbHhhpVWMwM/dwEvDkwcEb5v3xFw+5HEmHm9psKhyeTQtjPkmj1TefGpxNMHg+mU01Z3NWD7uBlaXlU4E97W/KzM3AZoCIuD8z14ynvNFq0rZIDWM2NWRbpAbqm09NzSZo3vboEK/WWA/3Aasi4oyIeB5wMbC94pokyWySVFfmkxrJmbMayMwDEXEFcAewBLg+M79QcVmSppzZJKmuzCc1lc1ZTWTmbcBtC/jI5lHVUoEmbYvUKGaTpLpaYD417e+5adujQmQecW63JEmSJGnMPOdMkiRJkmrA5mzCRMTaiNgZEbsiYmPV9fQTEddHxL6IeKA0tjQi7oyIB4vHE0uvXV1s286IuKCaqiUNwnySVEdmkyaJzdkEiYglwHXAm4DVwDsiYnW1VfW1BVjbNrYR2JGZq4AdxTLFtlwMnFV85gPFNkuqOfNJUh2ZTZo0NmeT5RxgV2Y+lJnPAjcC6yuuqafMvAvY3za8HthaPN8KXFgavzEzn8nMh4FdzG2zpPoznyTVkdmkiWJzNllOAR4rLe8uxibN8szcC1A8LivGm7J90jRqyt+v+SQ1S1P+ds2mKWFzNlmiw1iTLrfZ9O2Tmqzpf79N3z6pqZr+t9v07Zs6NmeTZTewsrR8KrCnoloW4/GIWAFQPO4rxpuyfdI0asrfr/kkNUtT/nbNpinhfc4mSEQcDfwd8PoXHnPMQy867jgeW7r4+4iv3H+g62vDWP+ofGPXvicz8+Sq65B0KJ+OPuHYM45ddkLV5VTOfJLqYRzZtHL/gVrvL5WZTfU3Gf8lCYDMPBARVwB3vOi44/j+rZfy/UNa97Xbnjhs+cqL5v5uh7X+Ubhn3aavVF2DpDmtfHrFU9/+xPP//J1Vl1M580mqh1Y2HbvshE+cvemnR/Z76ry/VGY21Z/N2YTJzNuA245ftTxbQXD3uk0AnH/rVQOts9yYtZoySVqozLzt+FXLObvqQiSppJVN0iSwOZtQL9u172BT1nL3uk0LbtBszCQNU6dsGvSLI0kaFrNJk8LmrCEGDZj5NGStBs7mTdJCufMjqY7MJtWVzVkDtAJmkJmzftrPRZOk+XDHR1IdmU2qO5uzCVcOmVEEjrNlkhbKnR9JdbPzzGWcP8ILgkjD4n3OJpg7QJLqZueZy6ouQZKkiWVzNqG67QC1n+wqSZIkaTJ4WGMPEXFuZn6631jVOl21EfrPrF277QmuvOjkjueVeTijpGEYxtXRvCiRJGla2Jz19kfAK+cxdpiIuB54M7AvM3+oGFsKbANOBx4B3p6ZXyteuxq4DHgO+JXMvGMhRZYvCFJeno9ujZk7Q1LzjDubYHGHX3fKp9aXSpKao4psWgy/1NYo2Zx1EBGvBn4UODkifq300gnAknmsYgvwH4A/K41tBHZk5vsiYmOx/K6IWA1cDJwFvAj4q4h4aWY+t9C6F7oTVA6SctDYmEmNtYUKsmlQrQxqZZJXj5UaawsTmE3SKNicdfY84HuZ++dzfGn8KeAn+304M++KiNPbhtcDry2ebwU+CbyrGL8xM58BHo6IXcA5wH8fuPoBlHeCDB2pmSYxm8AdIanpJjWbpFGwOesgMz8FfCoitmTmV4a02uWZubdY/96IaF3R4xTgb0rv212M9VS+0/0wr9roTpA0dYaaTTDYIdaS1Gbo2SRNApuz3v5DRGTb2NeB+4EPZea3h/A7osNY+++ce2PEBmADwGml8fKO0ChuRC1pKg01myRpSAbKpmNOPr7TW6Ta8VL6vT0EfAP4k+LnKeBx4KXF8kI8HhErAIrHfcX4bmBl6X2nAns6rSAzN2fmmsxc02l+y8voSxrAULNpf9ttPmzMJA1oqNn0PS98/kiLlYbF5qy3/zUzfyozby1+LgHOyczL6XPFxg62A5cWzy8FbimNXxwRx0TEGcAq4N5+K9t55rKp3umxEZWGZqjZBHMNWetHkgY09GySJoGHNfZ2ckSclpmPAkTEacBJxWvPdvtQRHyUuZNYT4qI3cC7gfcBN0XEZcCjwNsAMvMLEXET8EXgAHD5Qq44NK2HMp5/61VggyYtyDizSZLmy2ySDrE56+3XgXsi4v9j7hjnM4BfiojjmLtyUEeZ+Y4uL72+y/t/F/jdQYuctsZM0mDGnU2SNB9mk3SIhzV2EBE/ApCZtzE3Xf6rxc/LMvMTmfnNzPx3VdYI3vNHkiRJahKbs84+0HqSmc9k5t9m5meHdHXGoVi5/wBggyapXlbuP2AuSZI0IJszSdLQPLbUo+UlSRqU/xft7AciYnu3FzPzLeMsppPHlh49tTeMbl2psdONTiRVb1qzSZKkxbI56+wJ4A+rLkKdHbwAildrlCRJUoPYnHX2dGZ+quoiJEmSJE0Pzznr7JH2gYi4ZvxlTLe7123yZtOSJEmaGs6cdZCZb+0w/BbgmjGXMpXKDZn3cJMkSdK0cOZs/rz+xBi0z5Q5cyZJkqRp4czZ/M1UXUCTtZowZ8okSZI0rZw56yEifj8iToiI7wHujIgnI+KSquuSJEmS1DzOnPX2xsz8rYj4CWA38Dbgr4Ebqi2reZwxkyRJ0rRz5qy37ykefxz4aGbur7IYSZIkSc3lzFlv2yPiy8C3gF+KiJOBb1dckyRJkqQGcuasi4g4CrgVeDWwJjO/A/xPYH2lhUmSJElqJJuzLjLzu8AfZubXMvO5YuybmfnVikuTJEmS1EA2Z739ZUT8s4jwHmdD5L3LJEmSpCN5zllvvwYcBzwXEd9i7kbUmZknDLrCiHgEeBp4DjiQmWsiYimwDTgdeAR4e2Z+bXGl10+5Kbt73Sav0CjVyDRnk6R6M580TZw56yEzj8/MozLzezLzhGJ54Mas5HWZeXZmrimWNwI7MnMVsKNYbpxWM3b+rVfZmEn1NJXZJGkimE+aCs6c9VAczvjTwBmZ+a8jYiWwIjPvHfKvWg+8tni+Ffgk8K4h/47KtB/GWF62SZNqrdHZJGmimU9qJGfOevsAc1dr/Kli+RvAdYtcZzJ3LttsRGwoxpZn5l6A4nHZIn/Hgl277YmDP8PU6/wyGzOpVmqZTZKE+aQp4sxZb6/KzFdGxP8AyMyvRcTzFrnOczNzT0QsA+4s7qM2L0UgbQA45uTjF1nG4a686OShrq+bYTRkrQZyZtFrklRSy2ySJAbMJ7NJk8jmrLfvRMQS5r6xobgJ9XcXs8LM3FM87ouIm4FzgMcjYkVm7o2IFcC+Lp/dDGwGOH7V8lxMHePQbdZskMMahz2jJ+lw05RNw9Ypn/zySBqeQfNp2rOppZxRZlP9eVhjb/8euBlYFhG/C9wD/N6gK4uI4yLi+NZz4I3AA8B24NLibZcCtyym6Drrdf5ZSzlERnGopaTDmU2Hu3vdpoM//ZhP0miZT4MZ1ekqGj1nznrIzD+PiFng9cxdRv/CzPzSIla5HLi5uG3a0cBHMvP2iLgPuCkiLgMeBd62yNIrN997mZ1/61Udg8MwkcZqarKpl1ZutWb03bGRasF86uLabU9w5UUnzzunrrzoZLhhxEVp0WzOOoiIVzE3Df4S4PPAZZn5xcWuNzMfAl7RYfwfmGsAG+HudZuIa44cz2Ks9drsrkvc8ZFqYFqyqay1U9Nv9n6hDp6/6w6QNBTTmE8t8zkNpFdejet6Ahoum7POrgN+A7gLeAuwCbig0opqrPxtc7fGDOAzl1zCzJk3MLvrkoF/12FB486PpAHcM7uJmTOB2bkvi4aWSZK0CNdue4JX3nBo5+Yzl1xy2GvzZS5NNpuzzo7KzDuL5x+LiKsrraZmDjux9MwbDjVjs0c2Zu07PQvdCTJgJA1D60uk8sz9YphNkobh7nWbDn65DYc3ZAtlLjWDzVln3xcRb+22nJkfr6Cmsev1Lc3Mmf2nrfw2WtKozPdb5FZWlb84GjSbzCVJ81XOqG7ZcU/xpfbsticGnsU3l5rH5qyzTwHruiwn0JjmbKE7OKPQOpnVgJHUrpVR/U56b2XU7K5LOuaVOz2SRqlXPrVfuKOcVy2dMqqVQe3rNpuazeasg8z8uaprGJVrtz2x6PO+hqUcLgaNpLL2nZH5zuS3N2YeSi1plObz5fLMmTfALHPnuhbas6nciLU3ZebSdLE56yIiPpiZ74yI6zLz8qrrGUT7uWGtICgHwrBmxLrtABkokroZZOa+U9a0z5b1a8jMJUmL1amJ6qdXNnVqxMyq6WRz1kFEnAbcExHbgW0RcVpmPlp1XYvRKRAGbczqMOsmafL12vG4Z7bzvRJ75dZ8s8nDqCUtRquRah2NNB/lfOp0mLaZpBabs85eB5wGvBy4F1gC/FmlFfXQPkPW8sob5k6Cb7+/2GLMd+fHkJFU1qshar+S4iBa2dT+O9wBkrRY7fdDnG9WlXPpypnDX+t0GKMENmcdZebWiPhT4FXAezPzvVXX1G7l/gMHv7EpH8Nc1mrMhrnD0x4u4DHRkror79TcvcijqFtfNHW6GWunbAJzSVJ3nW5CD4cypv311v1ay1q5VP5cS7dcOuw9ZpTa2Jx192HmZs7eUx6MiLWZeXs1JR3ygv37+06llwNjvrrNjPX6ZsdgkVTWbYdnEOfNXHUwf86/dWirlTSl5pNPnd7Tarxmt801aOfNXFWMD7c+yeasg4j4FeBy4EvA2RFxZWbeUrz8e0DlzdnsikPP22fH2puyXrNnHqYoaZhetmvfUNYzu+uSQzP25o+kRXrZrn2L+uLo7nWb+Mwlc7l0HkfO3kvDYnPW2S8AM5n5jYg4HfjPEXF6Zl4LRKWVdVA+r2whyjs/7TxUUdKotB8y1OkG0WaPpCqUD01sZVSnQ6mlUbE562xJZn4DIDMfiYjXMtegvZiaNWftTdn5t17F3es2cf6tc4cCvfKGGw6+r7UD1Hr+yhtugIs6B447RpKGoddOTeu180pj8zlHQ5IGtZBGy6ZMVbA56+yrEXF2Zn4WoJhBezNwPXPnodVKe3i0lq+86OTDmq/zDr4+99zjpCUN284zl3H+pp+uugxJOozZpElhc9bZzwAHygOZeQD4mYj4UDUlHW5mLzx/5irOv7WZ0+69jgtvzQ7WagpTkiRJWiSbsw4yc3eP1z49zlp6aR2+OC1N2ULeI0mj0C9/mpTHkqTxszmbYE3ZCbDZkjTJmpLFkqTq2ZzVRESsBa4FlgB/mpnv6/X+nWcu4+yxVDZa823MWocyShqvhWZT09mISfVhPqmJbM5qICKWANcBbwB2A/dFxPbM/GK1lY1Wudn6zCW977d27bYnDrvapKTRm9ZsklR/5pOayuasHs4BdmXmQwARcSOwHmh0wLQu9z9fNmXS2E1lNkmaCOaTGsnmrB5OAR4rLe8GXtX+pojYAGwoFp+5Z92mB8ZQ20gVtzQ6CXhygI+/eJi1SDrC1GZTYdBsAvNJGrW++dTgbAL3nRrL5qweOl0VPo8YyNwMbAaIiPszc82oCxuHJm2L1DBmU0O2RWqgvvnU1GyC5m2PDjmq6gIEzH3bs7K0fCqwp6JaJKnFbJJUV+aTGsnmrB7uA1ZFxBkR8TzgYmB7xTVJktkkqa7MJzWShzXWQGYeiIgrgDuYuxzs9Zn5hT4f2zz6ysamSdsiNYbZ1KhtkRplgHxq2t9z07ZHhcg84vQBSZIkSdKYeVijJEmSJNWAzZkkSZIk1YDN2YSJiLURsTMidkXExqrr6Sciro+IfRHxQGlsaUTcGREPFo8nll67uti2nRFxQTVVSxqE+SSpjswmTRKbswkSEUuA64A3AauBd0TE6mqr6msLsLZtbCOwIzNXATuKZYptuRg4q/jM/9/evYdbUpV3Hv/+bC7iDaKgclPQ4AWNgwERRQ2ogy0zBjEi3iIYJ4qD95GowVHUREUm+owKCI8SMCqC8dbeQEQQUVEutkCLRCIoHRiJUYmogOA7f9Q69u7de+9zTp9Dnzqnv5/nOc+pXbWqatXatVbVW7Wq9nFt1QeuWwAAGOZJREFUmyX1nO2TpD6ybdJiY3C2uOwJXFVVP6qqW4GPAwcscJ4mqqrzgJ8PjT4AOKUNnwI8fWD8x6vqlqq6GriKbpsl9Z/tk6Q+sm3SomJwtrhsD1w78Hl1G7fY3Keqrgdo/+/dxi+V7ZM2Rkul/to+SUvLUqm7tk0bCYOzxSUjxi2l30JY6tsnLWVLvf4u9e2TlqqlXneX+vZtdAzOFpfVwI4Dn3cArlugvMzFT5NsC9D+39DGL5XtkzZGS6X+2j5JS8tSqbu2TRsJg7PF5UJglyQ7J9mM7gHQFQucp/WxAjikDR8CfHZg/LOTbJ5kZ2AX4DsLkD9Js2f7JKmPbJu0qGyy0BnQzFXVbUleBpwJLANOqqpVC5ytiZKcCuwDbJ1kNfBm4J3A6UleBPwEOAigqlYlOR34PnAbcHhV3b4gGZc0K7ZPkvrItkmLTarslipJkiRJC81ujZIkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDm8znwh673Xb1y1tumc9FShuVa+85r1VS2ujc9OsbFjoL0qK3+/ULnQNpcbsYzqyq5esz77yeCf7yllv4yP77z+cipY3KKw/eZqGzIC1q51/8noXOgrToXXTUQudAWtwCW6/vvHZrlCRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeMDiTJEmSpB4wOJMkSZKkHjA4kyRJkqQeSFXN38KSy4Gb522B0sZpa+BnC50JaRGzDklzYx2S5ubOVfXw9Zlxk3nOyM1Vtcc8L1PaqCS5yHokrT/rkDQ31iFpbpJctL7z2q1RkiRJknrA4EySJEmSemC+g7MT53l50sbIeiTNjXVImhvrkDQ3612H5vWFIJIkSZKk9WO3RkmSJEnqgTkFZ0kOSrIqye+TjH2rT5LlSa5MclWS189lndJSk+SeSc5K8sP2/4/GpLsmyWVJVs7lLUDSUjHdsSXJPklubHVmZZI3LUQ+pT5KclKSG9rPII2abv2RJkiyY5JzklzR4qFXjkiTJO9tx6lLk/zpdMud652zy4FnAOeNS5BkGXAs8FRgV+A5SXad43qlpeT1wNlVtQtwdvs8zr5VtZuvONbGbhbHlq+3OrNbVb11g2ZS6reTgeXTpLH+SOPdBvyvqnoosBdw+Ijj0FOBXdrfi4Hjp1vonIKzqrqiqq6cJtmewFVV9aOquhX4OHDAXNYrLTEHAKe04VOApy9gXqTFwmOLNAdVdR7w84XOh7RYVdX1VXVJG/4VcAWw/VCyA4APV+cCYKsk205a7oZ45mx74NqBz6tZN+PSxuw+VXU9dBUduPeYdAV8OcnFSV68wXIn9dNMjy2PSfK9JF9K8rANkzVpybD+SDOQZCfgkcC3hybNOg7aZAYr+wpw3xGTjqyqz043P5AR43xFpDYqk+rRLBazd1Vdl+TewFlJftCufEobo5kcWy4B7l9VNyXZH/gMXdcSSdOz/kgzkORuwCeBV1XVfw5PHjHLxDho2uCsqp488+yNtBrYceDzDsB1c1ymtKhMqkdJfppk26q6vt3qvmHMMq5r/29I8mm6bl0GZ9pYTXtsGTxIVtUXkxyXZOuq+tkGyqO0aFl/pOkl2ZQuMPtoVX1qRJJZx0EbolvjhcAuSXZOshnwbGDFBlivtFisAA5pw4cA69yRTnLXJHefGgb2o3shj7SxmvbYkuS+SdKG96Q75v3HBs+ptAhZf6TJWv34EHBFVb17TLIVwAvaWxv3Am6cepRlnGnvnE2TqQOB9wHbAF9IsrKqnpJkO+CDVbV/Vd2W5GXAmcAy4KSqWjWX9UpLzDuB05O8CPgJcBDAYD0C7gN8uh0nNwE+VlVnLFB+pQU37tiS5LA2/QPAM4GXJrkN+C3w7KqyW70EJDkV2AfYOslq4M3ApmD9kWZob+AvgcuSrGzj/ha4H/yhHn0R2B+4CvgN8MLpFhrrmSRJkiQtvA3RrVGSJEmSNA2DM0mSJEnqAYMzSZIkSeoBgzNJkiRJ6gGDM0mSJEnqAYMzzYskD0nyrSS3JHnthHSHJvn3JCuTfD/JX0+z3HOT7DH/OV5rHTslWdDfDEty0yzS7pzk20l+mOS09htPw2l2S7L//OZybH5eleQus0j/hCSXJLktyTOnSXtCkr2Hxm3etvuqVg47jZn3wW3/WZnkiiQnzjSPI5Z1nyQfS/KjJBe3ff3A9V1enyU5alIdHpH+De27uDLJU8akmdU+sr7WZ79Psrzl/aokr5+QbtskXx4x/qQkN8ykDUnymiQ/SHJZku8leXf7AdMlZbZtapK/T3LtpHYwyT5JHjs/OZw2P387y/QHJVmV5PfTHa+SnJFk+6Fxx7T94tIkn06y1Zh5T05ydWvTvpfkSbPM58nTtbl3pNm0LUnuleScJDclef+EdIu+bZGGGZxpvvwceAXwf2aQ9rSq2o3u91XenuQ+d0SGkszpd/x6vP6jgfdU1S7AL4AXjUizG93vaszYHPL7KmA2B8efAIcCH5tB2kcDFwyNexHwi6r6Y+A9dOUxynvpymm3qnoo3W8yzlq6H5f7DHBeVT2gqnan+8HjHUakXdB9bpx05r29T7IrXVk8DFgOHJdk2Yiks91HGLOc6cxqv2/rOBZ4KrAr8Jy2TaMsp/tNtWEnt2nTreswuh+P36uq/gR4FHADsMWYfPXOHbh/fw7Yc5o0+wCzCs7mkN9ZBWfA5cAzgPOmyc8WwD2r6t+GJp0FPLyqHgH8C/CGCYs5oh0/XwV8YJb5nHd3VNsC3Az8b2C6YG4ptC3SWgzONC+q6oaquhD43WzmAf4VuH+S45Nc1K4+vmVU+iT7tTsWlyT5RJK7jUhzbpK3J/ka8Mokuyf5WrvbcWaSbVu63duVx28Bh892e5Nck+ToJN9pf3/cxp/croafAxyd5IHtSunFSb6e5CEt3c5tWy5M8rZZrDfAE4F/bqNOAZ4+lGYz4K3Awe0K68FJ9kzyzSTfbf8f3NIe2sryc8CXk9wlyentCu5p6e5M7dHSrlP+SV4BbAec07Z5WlV1TVVdCvx+mm19KPAvVXX70KQD2nbTyuFJrVyGbQusHljvZW25y9qV6gvbdr6kjb9bkrPb9l2W5IA26xOBW9uPSU4t68dV9b4233AZ3jXd3ZQLW3kfMM1692n77T+nu3r+0THbM66cDk3y2bafXZnkzW38TunuGB4HXALsmOSIgfW/ZWAZR7Z5vwI8eKbrpvsuPl5Vt1TV1XQ/srnWSfaofWRcfW/16k1JzgcOSrJ/K5Pzk7w3yedbunXKeNR+P4P87wlcVVU/qqpbgY+3bRplOfCl4ZFVdR7dxanpHAm8tKp+2ea7tareWVX/2bbppiRvTfJt4DFJnt/alpXp7iAva+lGtoOt7N4ysP8+ZAZ5+oO2/n9o85+dZJs2/g5vU6vqgqq6fkLedgIOA17dyuPxSZ7W2qfvJvlK2kW+dHdnTkx3l/PDSbZJclbbrhOS/DjJ1i3tOmWc5J3AFm3cR2eY/yuq6soZJN0HOHfE/F+uqtvaxwsYceFnhG8B27ftGNe2JMn70/VS+QJw7xks9w8Wsm2pql9X1fl0Qdq4/C2VtkVaW1X559+8/QFHAa+dMP1Q4P1t+AF0V47vSXc1EWAZ3cHrEe3zucAewNZ0VyXv2sa/DnjTiOWfCxzXhjcFvgls0z4fDJzUhi8F/qwNHwNcPia/K8eMvwY4sg2/APh8Gz4Z+DywrH0+G9ilDT8a+GobXgG8oA0fDtzUhu8OrBzzt2srh6sG8rHjqLwPlnP7fA9gkzb8ZOCTA+lWD5T/a4ET2vDDgdumK/9WFlsPrOu0Mfl/wVAeTwaeOWFfeQ3wVyPGXw7sMPD5XwfXPzD+hcCNdCfUrwa2auNfDLyxDW8OXATsDGwC3KON35ou0AjdHeH3TLNPD5bh24Hnt+Gt6K6E33XCevdp+dyB7oLZt4DHjVjPYcBhY9Z/PXAvurswl7fvbCe6AHivlm4/4MS2TXei20+fAOwOXEZ39fkebbtf2+Y5Ysx3+d42/f1T29o+f2jUdzpiHxlX368B/qYN3xm4Fti5fT6VNfVsXBkfytr7/b5j8v/NNv2ZwAcH0v/l4PwD45cxpi1o03diTBsyUK9/MW56S1PAs9rwQ+nuJm3aPh9H185MVw9f3ob/5+B2DaxjO+CLE9b/vDb8Jta00+cyhzaV7oR8XJu21VAebppQPkcxcGwB/ghIG/4fwD8MpLsY2GJgH31DG17etnPrcWU8Kh/A18fk/8lD6c4F9piwDe8FnjjNfvA5BurU0LSTafWL7qLcx2pym/YMurtyy9p3/0tG18/etS1DeVinTi6VtsU//0b99bILjpa8g5M8DrgFeElV/TzJYUleTHeCvC1dIHLpwDx7tXHfSHdTYTO6k9hRTmv/H0wXYJzV5lkGXJ9kS7qTgq+1dP9E1/VgHdV1Hxnn1IH/7xkY/4mqur1d0X4s8ImsuRGyefu/N/AXA+s/uq3vV3TdJ0aaupo9nM0JeZyyJXBKkl1a+sHnXM6qqqkr/48D/m/Ly+VJpr6DGZd/Vc3kquJMPIUuwBo26q7SOmVQVf+Y5Ey6E7IDgJck+S90JxKPyJpnL7YEdqELsN6e5Al0Jx7bA+t0uU1yLF053VpVj2qjB8twP+DPs+bZijsD95uw3luB71TV6rb8lXQnP+cPbc+kLkxnVdV/tPk/1fL3GeDHVTXVLXS/9vfd9vlubf13Bz5dVb9p868YWOcxdCfa48zouxjhWRPq+1T9fQjwo+ruyEFXz148sC2jynjtjFSdw4T6NIv8Pxr49oTlTCeDy033bN7RdCd/z62qbwK3A59sSZ5Ed2J7YatvW9BdyJquHn6q/b+Y7sR8LVV1HeO7Zv2eNWX/kYFlwRza1OruKE36DtbXDsBp7c7dZsDVA9NWVNVv2/DjgANbXs5I8os2flwZr6OqHj9Ped6bCd30khxJd0Fs0h27Y5K8i+4u2F5t3Li25QnAqdX1PrguyVdHLbCnbcv6Wmxti7QOgzOtlySHA1Mv89i/HfSnTdf+n1ZVLxtIszPdAetRVfWLJCfTNYhrLYruIPGcGWTv1wPzrKqqxwzlaSvmp5GsMcNT678T8MsJAd46eUhyd7qrtKM8F7gC2CrJJtV1g9kBGFn2Q94GnFNVB7YuQueOyC+MPqBMjZ9R+Sc5jdHdV95dVR+eQV5J94D3VmP2q9V0dwxXp3umZEvg50n+HvhvsCaobvOfBJyU7gUFD2/b8vKqWuv5oSSHAtsAu1fV75JcQ7cfrmJNIE1VHd66RV00MPtwGf5FDXVzSncGOGq9+9BdqJhyO7Nvm4f3panPw/l6R1WdMLT+V42Yf2raEcDzRkw6r6pewZrvYsq0++MM6vtg/R27GEaX8aOHPu/L2hdOpvymqh47i/w/FThjQn7WzlyyI90dEIAPVNUHkvw6yc5VdXXbB85sXammXuhzc63pwhvglKp6w9Byn8bkeji1H63PPjRsVJs26zY1XRfq00ZNA/ap1s1zPbyPrk1Z0erQUSPyO5XnkVljRBmPTJh8nS7QGPbaqvrKTDKb5AHAtdV1cRs1/RDgvwNPqqpq4/4ReCRwXVVNHT+PoAucX0HXvXt3xrdp+zP3Y91CtS2ztkjbFmkdPnOm9VJVx1b3ooXdxgVms0h3D7pG88Z0zw2Muot1AbB31jzbdZckD5omm1cC2yR5TJtn0yQPaycDN7a7dzD6ADETBw/8X+cuUnXPklyd5KC2/rQ7NwDfoHuRwlrrr6pfDZTX8N/320H7HLouEwCHAJ9ty98zyVTw8yvWPpnYEph6CP3QCdt0PvCstrxdgT9p4yeV/1rrqqqDx+R/RoFZs2/bzlFW0G03dOXw1eocObWulsflaW/CS3Jfuq45/0b3UoeXDkx7UJK70pXRDS0w2xe4f1vHV4E7J3npQB4mPYB+JvDyFoyR5JED40etdz781yT3TPfCgafT7V+j8vVXWfOM0vZJ7k3XTe7AJFu0iwNPm5qhqo4Z811OnTytAJ6d7g2aO9NdLf9OW/6Hk0w9fza4j8ykvgP8AHhA1ryNc/CO7LgyHt4XzxmT/6kXS1wI7JLuGdDN6OrkCtb1JLouyjNSVdcOrGvqrsQ7gONbIDMVrA9fhJpyNvDM9v3Qvtv7s37t4EzdiTXtynMZunPbzLpNraorJ7RpEwOzJAcmeUf7OKlNO4TxBtu0/ei6Q8L4Mgb4XQbeollVjx+T/xkFZs3YAD/Jcrouqn8+dZeprfeFbT1r3e2sqt/T9XC4U7q7sOPalvPo6ueydHcY951FfqcsVNsy1hJrW6R1GJxpXiS5b5LVdM8JvTHJ6iT3mMm8VfU9uu4Qq+jucqzT+FfVv9MFFaem62p3AV3XhEnLvZXuZOPoJN+j6w8+1XC+EDg23cPrvx2ziKkuZuNsnu7h/VfSPdM0yvOAF7X1r2LNA8GvBA5PciHdScZsvA54TZKr6AKOD7Xx92PNtpwD7Jo1Dy+/C3hHkm/QdUUa5zi6k69L23ouBW6cpvxPBL6UGb4QJMmj2r5yEHBCklUjkk26U/Eh4F5t+18DjHtF8X7A5a3sz6R7y9n/Az4IfB+4JN3dtBPo7jJ8FNgjyUV039sPAFpA/HTgz9K9xvo7dFesXzdmvW+j6zZ6aVv+1Atfxq13RtJ1/T1szOTz6bqSraR7nvCi4QRV9WW6N2R+K8lldC9TuXtVXcKa5wQ/yfg7t+uoqlXA6W27zgAOH7j78wi651VgYB+ZSX1vy/4t3bNTZ6R7iP+ndM/mwfgyHt7vp8v/bcDL6PaPK4DT2zb9QbquxDe3iy3rSHIq3cWZB7d2b9TbUwGOB74CfLvVoW/QlcN3hxNW1feBN9K9YOZSuueGtl2fdnAor9sl+eKYyb8GHpbkYrqX4Lx1RL7m3KaOyde7Wptwl1aGR7VJDwSmyv1zdCf6K5M8nu5O2SfS3dX62YTFvwXYL8kldO3K9cCvxpVxm+dEun1rRi8EaUHkauAxwBfSdacetpzxbdr76U78z2rbN+1bGFu79HfA3zC+bfk08EO6576OB742all9bFtavq4B3g0c2vaLqbcdLom2RRpn6mFaSbPQDhp7VNWkk4INKskxwD9V9ybE9V3GMroH5G9O8kC6q8sPqjFdce4o7UTq0VU147d/bqzSdcfcowa6Ci+0dmHmQ1V10ByXc7equqldxT4W+GFVjepKdIdJ8ny6F9C8c0Oud0NLclNVrfMG3IWU5CPAq1tQur7L2By4vapua3f8jp+6u76htDx8o6ru0N/snG+2LdLC8JkzaYmoqiPmYTF3oXst8aZ0fe9fuqEDM4Cq+tMNvU7Nn3aXaU4nT81fp3sWZzO6K+InTJN+3lXVRzb0OtWpqufPw2LuB5ye7re4bmXNM9AbTFXdQveWQ83RUmpbpHG8cyZJkiRJPeAzZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AMGZ5IkSZLUAwZnkiRJktQDBmeSJEmS1AP/HzfRF4kt+jeoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x720 with 16 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from dataset import DatasetPrediction\n",
    "directory = os.getcwd()\n",
    "\n",
    "# Path to model, specify to load a model, otherwise let empty string\n",
    "model = r'weigths\\model30m.pth'\n",
    "model = os.path.join(directory,model)\n",
    "assert os.path.isfile(model), 'weight file not found'\n",
    "rain_dir = r'data\\pred_example\\rain'\n",
    "rain_dir = os.path.join(directory,rain_dir)\n",
    "assert os.path.isdir(rain_dir), 'rain file not found'\n",
    "wind_dir = r'data\\pred_example'\n",
    "wind_dir = os.path.join(directory,wind_dir)\n",
    "assert os.path.isdir(wind_dir), 'wind file not found'\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s')\n",
    "    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "    logging.info(f'Using device {device}')\n",
    "    \n",
    "    net = UNet(n_channels, n_classes, bilinear=True, n=8)\n",
    "    net.load_state_dict(\n",
    "        torch.load(model, map_location=device)\n",
    "    )\n",
    "    logging.info(f'Model loaded from {model}')\n",
    "    net.to(device=device)\n",
    "    net.eval()\n",
    "    \n",
    "    dataset = DatasetPrediction(rain_dir=rain_dir, \n",
    "                            U_dir=os.path.join(wind_dir,'U'), \n",
    "                            V_dir=os.path.join(wind_dir,'V'),\n",
    "                            thresholds=thresholds_in_cent_mm)\n",
    "    thresholds_normalized = np.log(np.array(thresholds_in_cent_mm)+1)/dataset.norm_factor\n",
    "    dataset_loader = DataLoader(dataset, batch_size=1, shuffle=False, num_workers=0, pin_memory=True)\n",
    "    batch = next(iter(dataset_loader))\n",
    "    imgs, true_imgs = batch['inputs'], batch['target']\n",
    "    imgs = imgs.to(device=device, dtype=torch.float32)\n",
    "    true_imgs = true_imgs.to(device=device, dtype=torch.float32)\n",
    "    persistance = utils.batch_to_mapped_persistance(imgs,thresholds_normalized)\n",
    "    # Network output\n",
    "    with torch.no_grad():\n",
    "        imgs_pred = net(imgs)\n",
    "        \n",
    "    utils.plot_comparison(imgs,true_imgs,imgs_pred,thresholds_normalized,temporal_length_inputs)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
